This notebook draws from previous notebooks on relative cumulative frequency analysis and PCA analysis of metabolic cage data.

Dataset

These datasets are from RMR data collected for mice after introduction to the metabolic cages (days 0-6, see image).

overview

A total of 24 mice were run over the length of 2 experiments.

Experiment 1: The first set of mice, run in February

Start of Met Cage run	Metabolic Cage #	Ear Tag	DOB	Gender	Genotype	Body Mass on BMR day (g)
2/4/2020	1	1049	10/08/19	M	WT	27.4
2/4/2020	2	1050	10/08/19	M	L-Bmal1-KO	26.6
2/4/2020	3	1053	10/07/19	M	WT	31.5
2/4/2020	4	1054	10/07/19	M	WT	30.5
2/4/2020	5	1055	10/07/19	M	L-Bmal1-KO	27.5
2/4/2020	6	1109	10/11/19	M	WT	28.3
2/4/2020	9	1110	10/11/19	M	WT	27.9
2/4/2020	10	1111	10/11/19	M	L-Bmal1-KO	26.7
2/4/2020	11	1112	10/11/19	M	WT	28.5
2/4/2020	12	1214	10/21/19	M	L-Bmal1-KO	25.3
2/4/2020	13	1215	10/21/19	M	WT	27.9
2/4/2020	14	1216	10/21/19	M	L-Bmal1-KO	27.6

Experiment 2: The second set of mice, run in March

Start of Met Cage run	Metabolic Cage #	Ear Tag	DOB	Gender	Genotype	Body Mass on BMR day (g)
3/10/2020	1	1555	11/18/19	M	WT	26.4
3/10/2020	2	1692	12/06/19	M	WT	25.6
3/10/2020	3	1693	12/06/19	M	WT	24.7
3/10/2020	4	1695	12/06/19	M	WT	25.6
3/10/2020	5	1699	12/06/19	M	WT	28.7
3/10/2020	6	1556	11/18/19	M	L-Bmal1-KO	26.2
3/10/2020	9	1557	11/18/19	M	L-Bmal1-KO	27.1
3/10/2020	10	1558	11/18/19	M	L-Bmal1-KO	26.9
3/10/2020	11	1559	11/18/19	M	L-Bmal1-KO	24.7
3/10/2020	12	1691	12/06/19	M	L-Bmal1-KO	28
3/10/2020	13	1694	12/06/19	M	L-Bmal1-KO	26.4
3/10/2020	14	1700	12/06/19	M	L-Bmal1-KO	27.7

Create data

Data generated using a python script.

SPF_vs_GF.py

#!/usr/bin/env python

"""
A CLI python script that uses the mousetrapper library to prepare WT vs KO mice data for import to R.
Must be run from same directory.
"""

import argparse
import numpy as np
import pandas as pd

import lib.mousetrapper as mt

__author__ = "Sumeed Yoyo Manzoor"
__email__ = "smanzoor@uchicago.edu"
__version__ = "0.0.1"
__experiment__ = "12 WT vs 12 L-Bmal1-KO mice"

datasets = "../datasets/"
data = "../data/"
outfolder = data + "R/"

experiment1_RMR = mt.Experiment.experiment(
    info="""First run with revised protocol, WT vs L-Bmal1-KO""",
    RT_data=datasets + "Bmal_WT_vs_KO_SPF_RMR_2020-02-04_rt.csv",
    macro13_data=datasets + "bmal_wt_vs_ko_spf_rmr_2020-02-04_macro_One-Click Macro V2.32.2-slice3-VOC.mac_1.csv", 
    t_ctrl={"SPF WT":[1,3,4,6,9,11,13]},
    t_exp={"SPF L-Bmal1-KO":[2,5,10,12,14]},
    groups={
        "SPF WT":[1,3,4,6,9,11,13],
        "SPF L-Bmal1-KO":[2,5,10,12,14]
    },
    meta=["WT","L-Bmal1-KO","WT","WT","L-Bmal1-KO","WT","WT","L-Bmal1-KO","WT","L-Bmal1-KO","WT","L-Bmal1-KO"],
    covar=[27.4,26.6,31.5,30.5,27.5,28.3,27.9,26.7,28.5,25.3,27.9,27.6]
    )

experiment2_RMR = mt.Experiment.experiment(
    info="""Second run with revised protocol, WT vs L-Bmal1-KO""",
    RT_data=datasets + "Bmal_SPF_KF_RMR_2020-03-10_14_23_rt.csv",
    macro13_data=datasets + "bmal_spf_kf_rmr_2020-03-10_14_23_macro_One-Click Macro V2.32.2-slice3-VOC.mac_1.csv",
    t_ctrl={"SPF WT":np.arange(1,6).tolist()},
    t_exp={"SPF L-Bmal1-KO":[6] + np.arange(9,15).tolist()},
    groups={
        "SPF WT":np.arange(1,6).tolist(),
        "SPF L-Bmal1-KO":[6] + np.arange(9,15).tolist()
    },
    meta=np.array(["WT","L-Bmal1-KO"]).repeat([5,7]).tolist(),
    covar=[26.4,25.6,24.7,25.6,28.7,26.2,27.1,26.9,24.7,28,26.4,27.7]
    )

experiment1_RMR.macro13_lightcycles(days=2)
experiment2_RMR.macro13_lightcycles(days=2)

experiment1_RMR.macro13_dark.to_csv(outfolder + "1_dark.csv", index=False)
experiment1_RMR.macro13_light.to_csv(outfolder + "1_light.csv", index=False)
pd.DataFrame({"covar":experiment1_RMR.covar}).to_csv(outfolder + "1_covar.csv", index=False)

experiment2_RMR.macro13_dark.to_csv(outfolder + "2_dark.csv", index=False)
experiment2_RMR.macro13_light.to_csv(outfolder + "2_light.csv", index=False)
pd.DataFrame({"covar":experiment2_RMR.covar}).to_csv(outfolder + "2_covar.csv", index=False)

shell("conda activate met_cages_development && cd py && python Bmal.py")

Import data

reqpkg <- c("magrittr","dplyr","ggplot2","ggpubr","DT","reshape2")
for (i in reqpkg) {
  print(i)
  print(packageVersion(i))
  suppressWarnings(suppressMessages(library(i, quietly=T, verbose=T, warn.conflicts=F,character.only=T)))
}

## [1] "magrittr"
## [1] '1.5'
## [1] "dplyr"
## [1] '0.8.5'
## [1] "ggplot2"
## [1] '3.3.0'
## [1] "ggpubr"
## [1] '0.2.5'
## [1] "DT"
## [1] '0.13'
## [1] "reshape2"
## [1] '1.4.3'

theme_set(theme_pubr())
select <- dplyr::select

path <- "data/R/"

exp1 <- read.csv2(paste0(path, "1_total.csv"), sep = ",") %>% 
  set_colnames(paste0("exp1_",colnames(.)))
exp1_dark <- read.csv2(paste0(path, "1_dark.csv"), sep = ",") %>% 
  set_colnames(paste0("exp1_",colnames(.))) %>% slice(-c(241, 242, 483, 484))
exp1_light <- read.csv2(paste0(path, "1_light.csv"), sep = ",") %>% 
  set_colnames(paste0("exp1_",colnames(.)))
exp1_covar <- read.csv2(paste0(path, "1_covar.csv"), sep = ",") %>% 
  mutate_all(function(x) as.numeric(as.character(x)))

exp2 <- read.csv2(paste0(path, "2_total.csv"), sep = ",") %>% 
  set_colnames(paste0("exp2_",colnames(.)))
exp2_dark <- read.csv2(paste0(path, "2_dark.csv"), sep = ",") %>%
  set_colnames(paste0("exp2_",colnames(.))) %>% slice(-c(241, 482))
exp2_light <- read.csv2(paste0(path, "2_light.csv"), sep = ",") %>% 
  set_colnames(paste0("exp2_",colnames(.))) %>% slice(-c(241, 482))
exp2_covar <- read.csv2(paste0(path, "2_covar.csv"), sep = ",") %>% 
  mutate_all(function(x) as.numeric(as.character(x)))

covar <- c(exp1_covar[c(1,3,4,6,7,9,11),], exp2_covar[1:5,], exp1_covar[c(2,5,8,10,12),], exp2_covar[6:12,])

total2 <- exp1 %>% 
  bind_cols(exp2) %>% 
  select(-contains(c("Enviro", "Time"))) %>% 
  mutate_all(function(x) as.numeric(as.character(x)))

light <- exp1_light %>% 
  bind_cols(exp2_light) %>% 
  select(-contains(c("Enviro", "Time"))) %>% 
  mutate_all(function(x) as.numeric(as.character(x)))

dark <- exp1_dark %>% 
  bind_cols(exp2_dark) %>% 
  select(-contains(c("Enviro", "Time"))) %>% 
  mutate_all(function(x) as.numeric(as.character(x)))

total <- bind_rows(light, dark)

Datatables

Light

Hide

Select a tab above to view

VO2

light %>% select(contains("VO2")) %>% datatable(filter="bottom")

RQ

light %>% select(contains("RQ")) %>% datatable(filter="bottom")

BodyMass

light %>% select(contains("BodyMass")) %>% datatable(filter="bottom")

Food

light %>% select(contains("Food")) %>% datatable(filter="bottom")

Water

light %>% select(contains("Water")) %>% datatable(filter="bottom")

PedMeters

light %>% select(contains("PedMeters")) %>% datatable(filter="bottom")

Dark

Hide

Select a tab above to view

VO2

dark %>% select(contains("VO2")) %>% datatable(filter="bottom")

RQ

dark %>% select(contains("RQ")) %>% datatable(filter="bottom")

BodyMass

dark %>% select(contains("BodyMass")) %>% datatable(filter="bottom")

Food

dark %>% select(contains("Food")) %>% datatable(filter="bottom")

Water

dark %>% select(contains("Water")) %>% datatable(filter="bottom")

PedMeters

dark %>% select(contains("PedMeters")) %>% datatable(filter="bottom")

Graphs

select_WT <- function(df) {
  d1 <- df %>% select(contains("exp1"))
  for (name in colnames(d1)) {
    if (as.numeric(strsplit(name, "_")[[1]][4]) %in% c(2,5,10,12,14)) {
      d1 %<>% select(-all_of(name))
    }
  }
  
  d2 <- df %>% select(contains("exp2"))
  for (name in colnames(d2)) {
    if (as.numeric(strsplit(name, "_")[[1]][4]) > 5) {
      d2 %<>% select(-all_of(name))
    }
  }
  
  d <- bind_cols(d1, d2)
  return(d)
}

select_KO <- function(df) {
  d1 <- df %>% select(contains("exp1"))
  for (name in colnames(d1)) {
    if (as.numeric(strsplit(name, "_")[[1]][4]) %in% c(1,3,4,6,9,11,13)) {
      d1 %<>% select(-all_of(name))
    }
  }
  
  d2 <- df %>% select(contains("exp2"))
  for (name in colnames(d2)) {
    if (as.numeric(strsplit(name, "_")[[1]][4]) <= 5) {
      d2 %<>% select(-all_of(name))
    }
  }
  
  d <- bind_cols(d1, d2)
  return(d)
}

graph_groups <- function(data, channel, title="", type="m") {
  df <- data %>% select(contains(channel))
  df.observation <- 1:nrow(df)
  c <- ncol(df)
  WT <- df %>% select_WT() %>% set_colnames(paste0("WT_", colnames(.))) %>% melt() %>% select(-variable)
  KO <- df %>% select_KO() %>% set_colnames(paste0("KO_", colnames(.))) %>% melt() %>% select(-variable)
  df <- bind_cols(WT, KO) %>% set_colnames(c("WT","KO")) %>% melt()
  df$observation <- rep(df.observation, c)
  
  if (type=="m") {
    df %>% group_by(variable, observation) %>% 
      summarise(value=mean(value)) %>% 
      ggplot(aes(x=observation, y=value, color=variable)) +
        geom_rect(aes(xmin=241, xmax=481, ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
        geom_rect(aes(xmin=741, xmax=961, ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
        geom_line(size=1) +
        ggtitle(title) +
        ylab(channel)
  } else if (type=="s") {
    ggplot(df, aes(x=observation, y=value, color=variable)) + 
      ggtitle(title) + 
      geom_rect(aes(xmin=241, xmax=481, ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
      geom_rect(aes(xmin=741, xmax=961, ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
      stat_summary(fun.data = "mean_cl_boot", geom = "errorbar", alpha = 0.7) +
      ylab(channel)
  }
}

graph_mice <- function(data, channel, title="") {
  df <- data %>% select(contains(channel))
  df.observation <- 1:nrow(df)
  c <- ncol(df)
  WT <- df %>% select_WT() %>% set_colnames(paste0("WT_", colnames(.))) %>% melt()
  KO <- df %>% select_KO() %>% set_colnames(paste0("KO_", colnames(.))) %>% melt()
  df <- bind_rows(WT, KO) %>% set_colnames(c("mouse", "channel"))
  df$observation <- rep(df.observation, c)
  ggplot(df, aes(x=observation, y=channel, color=mouse)) + 
    geom_rect(aes(xmin=241, xmax=481, ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
    geom_rect(aes(xmin=741, xmax=961, ymin=-Inf, ymax=Inf), fill="gainsboro", inherit.aes = FALSE) +
    geom_line() +
    ylab(channel)
}

Total cycle

Hide

Select a tab above to view

VO2

graph_groups(total2, "VO2", title="VO2 by group", type = "m")

## No id variables; using all as measure variables
## No id variables; using all as measure variables
## No id variables; using all as measure variables

graph_groups(total2, "VO2", title="VO2 by group", type = "s")

## No id variables; using all as measure variables
## No id variables; using all as measure variables
## No id variables; using all as measure variables

graph_mice(total2, "VO2", title="VO2 by mouse")

## No id variables; using all as measure variables
## No id variables; using all as measure variables

## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character

## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector

RQ

graph_groups(total2, "RQ", title="RQ by group", type = "m")

## No id variables; using all as measure variables
## No id variables; using all as measure variables
## No id variables; using all as measure variables

graph_groups(total2, "RQ", title="RQ by group", type = "s")

## No id variables; using all as measure variables
## No id variables; using all as measure variables
## No id variables; using all as measure variables

graph_mice(total2, "RQ", title="RQ by mouse")

## No id variables; using all as measure variables
## No id variables; using all as measure variables

## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character

## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector

BodyMass

graph_groups(total2, "BodyMass", title="BodyMass by group", type = "m")

## No id variables; using all as measure variables
## No id variables; using all as measure variables
## No id variables; using all as measure variables

graph_groups(total2, "BodyMass", title="BodyMass by group", type = "s")

## No id variables; using all as measure variables
## No id variables; using all as measure variables
## No id variables; using all as measure variables

graph_mice(total2, "BodyMass", title="BodyMass by mouse")

## No id variables; using all as measure variables
## No id variables; using all as measure variables

## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character

## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector

Food

graph_groups(total2, "Food", title="Food by group", type = "m")

## No id variables; using all as measure variables
## No id variables; using all as measure variables
## No id variables; using all as measure variables

graph_groups(total2, "Food", title="Food by group", type = "s")

## No id variables; using all as measure variables
## No id variables; using all as measure variables
## No id variables; using all as measure variables

graph_mice(total2, "Food", title="Food by mouse")

## No id variables; using all as measure variables
## No id variables; using all as measure variables

## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character

## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector

Water

graph_groups(total2, "Water", title="Water by group", type = "m")

## No id variables; using all as measure variables
## No id variables; using all as measure variables
## No id variables; using all as measure variables

graph_groups(total2, "Water", title="Water by group", type = "s")

## No id variables; using all as measure variables
## No id variables; using all as measure variables
## No id variables; using all as measure variables

graph_mice(total2, "Water", title="Water by mouse")

## No id variables; using all as measure variables
## No id variables; using all as measure variables

## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character

## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector

PedMeters

graph_groups(total2, "PedMeters", title="PedMeters by group", type = "m")

## No id variables; using all as measure variables
## No id variables; using all as measure variables
## No id variables; using all as measure variables

graph_groups(total2, "PedMeters", title="PedMeters by group", type = "s")

## No id variables; using all as measure variables
## No id variables; using all as measure variables
## No id variables; using all as measure variables

graph_mice(total2, "PedMeters", title="PedMeters by mouse")

## No id variables; using all as measure variables
## No id variables; using all as measure variables

## Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character

## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector

## Warning in bind_rows_(x, .id): binding character and factor vector, coercing
## into character vector

RCF

Libraries

reqpkg <- c("docstring","tibble","scales","purrr","nplr", "car","compute.es","effects","multcomp","pastecs","WRS2")
for (i in reqpkg) {
  print(i)
  print(packageVersion(i))
  suppressWarnings(suppressMessages(library(i, quietly=T, verbose=T, warn.conflicts=F,character.only=T)))
}

## [1] "docstring"
## [1] '1.0.0'
## [1] "tibble"
## [1] '2.1.3'
## [1] "scales"
## [1] '1.1.0'
## [1] "purrr"
## [1] '0.3.3'
## [1] "nplr"
## [1] '0.1.7'
## [1] "car"
## [1] '3.0.7'
## [1] "compute.es"
## [1] '0.2.4'
## [1] "effects"
## [1] '4.1.4'
## [1] "multcomp"
## [1] '1.4.12'
## [1] "pastecs"
## [1] '1.3.21'
## [1] "WRS2"
## [1] '1.0.0'

Functions

signif.floor <- function(x, n){
  pow <- floor( log10( abs(x) ) ) + 1 - n
  y <- floor(x / 10 ^ pow) * 10^pow
  # handle the x = 0 case
  y[x==0] <- 0
  y
}

signif.ceiling <- function(x, n){
  pow <- floor( log10( abs(x) ) ) + 1 - n
  y <- ceiling(x / 10 ^ pow) * 10^pow
  # handle the x = 0 case
  y[x==0] <- 0
  y
}

find_range <- function(data, channel, segments=10) {
  data <- data %>% dplyr::select(contains(channel)) %>% melt()
  r <- range(data$value)
  cut <- max(data$value)-min(data$value)
  out <- list("range"=r, "cut"=cut/segments)
  return(out)
}

rcf <- function(data, channel, id, min, max, b=0.1) {
  #' Relative cumulative frequency function
  #'
  #' Takes a vector and returns a dataframe with relative cumulative frequency per break unit
  #'
  #' @param data numeric vector
  #' @param channel character. Channel to label the rcf output.
  #' @param b double. Number to cut vector.
  #' Before running, check for valid break values to cut the data. One way is to divide range of data by 10 `(max(data)-min(data))/10`.
  
  breaks <- seq(min, max, by=b)
  # breaks <- seq(signif.floor(min(data), 2), signif.ceiling(max(data), 2), by=b)
  duration.cut <- cut(data, breaks, right = F)
  duration.freq <- table(duration.cut)
  duration.cumfreq = cumsum(duration.freq) %>% prepend(0)
  duration.cumrelfreq = duration.cumfreq/length(data)
  out <- duration.cumrelfreq %>% as.data.frame()
  out <- cbind(breaks, out) %>% set_colnames(c("breaks", channel)) %>% set_rownames(NULL)
  out$id <- rep(id, nrow(out))
  return(out)
}

EC50 <- function(x) {
  out <- lapply(x, function(y) {
    nplr(y$breaks, y[, 2], silent = T, useLog = FALSE)
  }) %>% 
    sapply(function(z) {
      getEstimates(z, 0.5)$x
    })
  out %<>% set_rownames(NULL)
  return(out)
}

lm_eqn <- function(df, y, x){
  m <- lm(y ~ x, df);
  eq <- substitute(italic(y) == a + b %.% italic(x)*","~~italic(r)^2~"="~r2, 
                   list(a = format(unname(coef(m)[1]), digits = 2),
                        b = format(unname(coef(m)[2]), digits = 2),
                        r2 = format(summary(m)$r.squared, digits = 3)))
  as.character(as.expression(eq));
}

lm_eqn_str <- function(df, y, x){
  m <- lm(y ~ x, df);
  eq <- sprintf("y = %s + %sx, r^2 = %s", format(unname(coef(m)[1]), digits = 2),format(unname(coef(m)[2]), digits = 2),format(summary(m)$r.squared, digits = 3))
  return(eq)
}

remove_outliers <- function(x, na.rm = TRUE, ...) {
  qnt <- quantile(x, probs=c(.25, .75), na.rm = na.rm, ...)
  H <- 1.5 * IQR(x, na.rm = na.rm)
  y <- x
  y[x < (qnt[1] - H)] <- NA
  y[x > (qnt[2] + H)] <- NA
  y
}

ANCOVA <- function(df) {
  # Takes dataframe with column 1 as independent variable (e.g., WT), column 2 as dependent variable, and column 3 as covariate
  
  cat("linear fit of dependent variable with covariates\n")
  p1 <- ggplot(df, aes_string(x=names(df)[2], y=names(df)[3])) +
    geom_smooth(method = "lm", se = FALSE) +
    geom_point() +
    # geom_text(x = mean(na.omit(df[[2]])), y = mean(na.omit(df[[3]])), label = lm_eqn(df, df[[3]], df[[2]]), parse = TRUE) +
    geom_text(x = mean(na.omit(df[[2]])), y = (mean(na.omit(df[[3]])) + IQR(df[[3]], na.rm = T)), label = lm_eqn(df, df[[3]], df[[2]]), parse = TRUE) +
    ggtitle("Linear model - DV and Covar") +
    xlab("DV") + ylab("Covar")
  print(p1)
  
  m <- lm(df[[3]] ~ df[[2]], df)
  
  summary(m) %>% print()
  lm_eqn_str(df, df[[3]], df[[2]]) %>% cat()
  
  plot1 <- ggplot(df, aes_string(x=names(df)[1], y=names(df)[2])) +
    geom_boxplot() +
    # ggtitle("DV vs IV") +
    xlab("IV") + ylab("DV")
  plot2 <- ggplot(df, aes_string(x=names(df)[1], y=names(df)[3])) +
    geom_boxplot() +
    # ggtitle("Covar vs IV") +
    xlab("IV") + ylab("Covar")
  p2 <- ggarrange(plot1, plot2,
                  labels = c("DV vs IV", "Covar vs IV"),
                  ncol = 2)
  # grid.arrange(plot1, plot2, ncol=2)
  print(p2)
  
  cat("\n_____________________\n")
  cat("Some stats on the variables\n")
  cat("DV - IV:\n")
  by(df[[2]], df[[1]], stat.desc) %>% print()
  cat("Covar - IV:\n")
  by(df[[3]], df[[1]], stat.desc) %>% print()
  cat("levene's test\n")
  leveneTest(df[[2]], df[[1]], center = median) %>% print()
  
  cat("_____________________\n")
  cat("ANOVA of DV with IV\n")
  model <- aov(df[[2]] ~ df[[1]], df)
  summary(model) %>% print()
  
  cat("_____________________\n")
  cat("ANOVA of covariate with IV, independent of DV\n")  
  model <- aov(df[[3]] ~ df[[1]], df)
  summary(model) %>% print()
  
  cat("_____________________\n")
  cat("ANOVA of DV with IV, normalized to covariate\n")
  c <- df[[2]]/df[[3]]
  model <- aov(c ~ df[[1]], df)
  summary(model) %>% print()
  
  cat("_____________________\n")
  cat("ANCOVA\n")
  model <- aov(df[[2]] ~ df[[1]] + df[[3]], data = df)
  summary(model) %>% print()
  out <- Anova(model, type= "III")
  print(out)
  return(out)
}

check_ANCOVA <- function(res) {
  if (res$`Pr(>F)`[2] < 0.05 && res$`Pr(>F)`[3] < 0.05) {
    print(paste("Varies significantly with DV and Covar,"), res$`Pr(>F)`[2], res$`Pr(>F)`[3], sep=" ")
  } else if (res$`Pr(>F)`[2] < 0.05 && res$`Pr(>F)`[3] > 0.05) {
    print(paste("Varies significantly with DV,", "p =", res$`Pr(>F)`[2], sep =" "))
  } else if (res$`Pr(>F)`[2] > 0.05 && res$`Pr(>F)`[3] < 0.05) {
    print(paste("Varies significantly with Covar,", "p =", res$`Pr(>F)`[3], sep=" "))
  }
}

VO2

Hide

Select a tab above to view

Total cycle

total.VO2 <- total %>% select(contains("VO2"))
total.VO2.WT <- total.VO2 %>% select_WT()
total.VO2.KO <- total.VO2 %>% select_KO()

rcf.WT <- lapply(total.VO2.WT, function(x) rcf(x, "VO2", "WT", 0.9, 3.2, 0.1))
rcf.KO <- lapply(total.VO2.KO, function(x) rcf(x, "VO2", "KO", 0.9, 3.2, 0.1))

data <- bind_rows(bind_rows(rcf.WT), bind_rows(rcf.KO)) %>% split(.$id)

Models <- lapply(data, function(x){
  nplr(x$breaks, x$VO2, silent = TRUE, useLog = F)
})

overlay(Models, xlab = expression(VO[2]~(mL/min)), ylab = "relative cumulative frequency", main="RCF of WT vs KO mice VO2 in complete cycle", cex.main=1.5, lwd = 3, Cols = hue_pal()(2))

print(paste0("WT EC50: ", getEstimates(Models$WT, 0.5)$x))

## [1] "WT EC50: 1.53037412338406"

print(paste0("KO EC50: ", getEstimates(Models$KO, 0.5)$x))

## [1] "KO EC50: 1.64418607579127"

t <- list(rcf.WT, rcf.KO) %>% sapply(EC50) %>% set_rownames(NULL) %>% set_colnames(c("WT", "KO"))
var.test(t[,1], t[,2])

## 
##  F test to compare two variances
## 
## data:  t[, 1] and t[, 2]
## F = 1.1019, num df = 11, denom df = 11, p-value = 0.8751
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.3172068 3.8276027
## sample estimates:
## ratio of variances 
##           1.101881

t.test(t[,1], t[,2])

## 
##  Welch Two Sample t-test
## 
## data:  t[, 1] and t[, 2]
## t = -1.5795, df = 21.948, p-value = 0.1285
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.24455700  0.03310949
## sample estimates:
## mean of x mean of y 
##  1.554837  1.660560

df <- data.frame(meta=c(rep("WT", nrow(t)), rep("KO", nrow(t))), ch=c(t[,1],t[,2]), covar=covar)
res <- ANCOVA(df)

## linear fit of dependent variable with covariates

## 
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.5032 -0.9561 -0.1400  0.6464  4.4487 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   25.605      3.295   7.772 9.52e-08 ***
## df[[2]]        1.016      2.039   0.498    0.623    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.654 on 22 degrees of freedom
## Multiple R-squared:  0.01116,    Adjusted R-squared:  -0.03379 
## F-statistic: 0.2482 on 1 and 22 DF,  p-value: 0.6233
## 
## y = 26 + 1x, r^2 = 0.0112

## 
## _____________________
## Some stats on the variables
## DV - IV:
## df[[1]]: KO
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  12.00000000   0.00000000   0.00000000   1.37522815   2.02098097   0.64575282 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##  19.92672416   1.65187506   1.66056035   0.04616875   0.10161673   0.02557864 
##      std.dev     coef.var 
##   0.15993323   0.09631281 
## ------------------------------------------------------------ 
## df[[1]]: WT
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  12.00000000   0.00000000   0.00000000   1.31808933   1.89266308   0.57457375 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##  18.65803908   1.56089558   1.55483659   0.04846357   0.10666760   0.02818461 
##      std.dev     coef.var 
##   0.16788274   0.10797452 
## Covar - IV:
## df[[1]]: KO
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  12.00000000   0.00000000   0.00000000  24.70000000  28.00000000   3.30000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 320.70000000  26.80000000  26.72500000   0.28394542   0.62495965   0.96750000 
##      std.dev     coef.var 
##   0.98361578   0.03680508 
## ------------------------------------------------------------ 
## df[[1]]: WT
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  12.00000000   0.00000000   0.00000000  24.70000000  31.50000000   6.80000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 333.00000000  27.90000000  27.75000000   0.57689083   1.26972816   3.99363636 
##      std.dev     coef.var 
##   1.99840846   0.07201472 
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  0.7281 0.4027
##       22               
## _____________________
## ANOVA of DV with IV
##             Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]]      1 0.0671 0.06707   2.495  0.128
## Residuals   22 0.5914 0.02688               
## _____________________
## ANOVA of covariate with IV, independent of DV
##             Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]]      1   6.30   6.304   2.541  0.125
## Residuals   22  54.57   2.481               
## _____________________
## ANOVA of DV with IV, normalized to covariate
##             Df    Sum Sq   Mean Sq F value Pr(>F)  
## df[[1]]      1 0.0002107 2.107e-04   6.008 0.0227 *
## Residuals   22 0.0007716 3.507e-05                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## _____________________
## ANCOVA
##             Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]]      1 0.0671 0.06707   2.517  0.128
## df[[3]]      1 0.0319 0.03188   1.196  0.286
## Residuals   21 0.5595 0.02664               
## Anova Table (Type III tests)
## 
## Response: df[[2]]
##              Sum Sq Df F value  Pr(>F)  
## (Intercept) 0.07817  1  2.9337 0.10147  
## df[[1]]     0.09160  1  3.4378 0.07782 .
## df[[3]]     0.03188  1  1.1964 0.28643  
## Residuals   0.55952 21                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

check_ANCOVA(res)

Light cycle

light.VO2 <- light %>% select(contains("VO2"))
light.VO2.WT <- light.VO2 %>% select_WT()
light.VO2.KO <- light.VO2 %>% select_KO()

rcf.WT <- lapply(light.VO2.WT, function(x) rcf(x, "VO2", "WT", 0.9, 3.2, 0.1))
rcf.KO <- lapply(light.VO2.KO, function(x) rcf(x, "VO2", "KO", 0.9, 3.2, 0.1))

data <- bind_rows(bind_rows(rcf.WT), bind_rows(rcf.KO)) %>% split(.$id)

Models <- lapply(data, function(x){
  nplr(x$breaks, x$VO2, silent = TRUE, useLog = F)
})

overlay(Models, xlab = expression(VO[2]~(mL/min)), ylab = "relative cumulative frequency", main="RCF of WT vs KO mice VO2 in light cycle", cex.main=1.5, lwd = 3, Cols = hue_pal()(2))

print(paste0("WT EC50: ", getEstimates(Models$WT, 0.5)$x))

## [1] "WT EC50: 1.38025535607356"

print(paste0("KO EC50: ", getEstimates(Models$KO, 0.5)$x))

## [1] "KO EC50: 1.44911832907739"

t <- list(rcf.WT, rcf.KO) %>% sapply(EC50) %>% set_rownames(NULL) %>% set_colnames(c("WT", "KO"))
var.test(t[,1], t[,2])

## 
##  F test to compare two variances
## 
## data:  t[, 1] and t[, 2]
## F = 0.70246, num df = 11, denom df = 11, p-value = 0.568
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.2022238 2.4401508
## sample estimates:
## ratio of variances 
##          0.7024646

t.test(t[,1], t[,2])

## 
##  Welch Two Sample t-test
## 
## data:  t[, 1] and t[, 2]
## t = -1.2306, df = 21.348, p-value = 0.2319
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.22122497  0.05663625
## sample estimates:
## mean of x mean of y 
##  1.381931  1.464225

df <- data.frame(meta=c(rep("WT", nrow(t)), rep("KO", nrow(t))), ch=c(t[,1],t[,2]), covar=covar)
res <- ANCOVA(df)

## linear fit of dependent variable with covariates

## 
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.4395 -1.0950 -0.1280  0.6293  4.5810 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   23.583      2.895   8.145 4.37e-08 ***
## df[[2]]        2.568      2.021   1.271    0.217    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.606 on 22 degrees of freedom
## Multiple R-squared:  0.06837,    Adjusted R-squared:  0.02602 
## F-statistic: 1.614 on 1 and 22 DF,  p-value: 0.2171
## 
## y = 24 + 2.6x, r^2 = 0.0684

## 
## _____________________
## Some stats on the variables
## DV - IV:
## df[[1]]: KO
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  12.00000000   0.00000000   0.00000000   1.19716087   1.91023222   0.71307134 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##  17.57070207   1.43327793   1.46422517   0.05125155   0.11280390   0.03152066 
##      std.dev     coef.var 
##   0.17754058   0.12125224 
## ------------------------------------------------------------ 
## df[[1]]: WT
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  12.00000000   0.00000000   0.00000000   1.17764983   1.62562865   0.44797882 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##  16.58316975   1.37143501   1.38193081   0.04295555   0.09454452   0.02214215 
##      std.dev     coef.var 
##   0.14880238   0.10767715 
## Covar - IV:
## df[[1]]: KO
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  12.00000000   0.00000000   0.00000000  24.70000000  28.00000000   3.30000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 320.70000000  26.80000000  26.72500000   0.28394542   0.62495965   0.96750000 
##      std.dev     coef.var 
##   0.98361578   0.03680508 
## ------------------------------------------------------------ 
## df[[1]]: WT
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  12.00000000   0.00000000   0.00000000  24.70000000  31.50000000   6.80000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 333.00000000  27.90000000  27.75000000   0.57689083   1.26972816   3.99363636 
##      std.dev     coef.var 
##   1.99840846   0.07201472 
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1       0 0.9945
##       22               
## _____________________
## ANOVA of DV with IV
##             Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]]      1 0.0406 0.04063   1.514  0.231
## Residuals   22 0.5903 0.02683               
## _____________________
## ANOVA of covariate with IV, independent of DV
##             Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]]      1   6.30   6.304   2.541  0.125
## Residuals   22  54.57   2.481               
## _____________________
## ANOVA of DV with IV, normalized to covariate
##             Df    Sum Sq   Mean Sq F value Pr(>F)  
## df[[1]]      1 0.0001430 1.430e-04   4.812 0.0391 *
## Residuals   22 0.0006536 2.971e-05                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## _____________________
## ANCOVA
##             Df Sum Sq Mean Sq F value Pr(>F)  
## df[[1]]      1 0.0406 0.04063   1.682 0.2088  
## df[[3]]      1 0.0829 0.08287   3.429 0.0782 .
## Residuals   21 0.5074 0.02416                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
## 
## Response: df[[2]]
##              Sum Sq Df F value  Pr(>F)  
## (Intercept) 0.01357  1  0.5618 0.46187  
## df[[1]]     0.08037  1  3.3260 0.08246 .
## df[[3]]     0.08287  1  3.4294 0.07816 .
## Residuals   0.50743 21                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

check_ANCOVA(res)

Dark cycle

dark.VO2 <- dark %>% select(contains("VO2"))
dark.VO2.WT <- dark.VO2 %>% select_WT()
dark.VO2.KO <- dark.VO2 %>% select_KO()

rcf.WT <- lapply(dark.VO2.WT, function(x) rcf(x, "VO2", "WT", 0.9, 3.2, 0.1))
rcf.KO <- lapply(dark.VO2.KO, function(x) rcf(x, "VO2", "KO", 0.9, 3.2, 0.1))

data <- bind_rows(bind_rows(rcf.WT), bind_rows(rcf.KO)) %>% split(.$id)

Models <- lapply(data, function(x){
  nplr(x$breaks, x$VO2, silent = TRUE, useLog = F)
})

overlay(Models, xlab = expression(VO[2]~(mL/min)), ylab = "relative cumulative frequency", main="RCF of WT vs KO mice VO2 in dark cycle", cex.main=1.5, lwd = 3, Cols = hue_pal()(2))

print(paste0("WT EC50: ", getEstimates(Models$WT, 0.5)$x))

## [1] "WT EC50: 1.76003338341844"

print(paste0("KO EC50: ", getEstimates(Models$KO, 0.5)$x))

## [1] "KO EC50: 1.8830318193838"

t <- list(rcf.WT, rcf.KO) %>% sapply(EC50) %>% set_rownames(NULL) %>% set_colnames(c("WT", "KO"))
var.test(t[,1], t[,2])

## 
##  F test to compare two variances
## 
## data:  t[, 1] and t[, 2]
## F = 1.4877, num df = 11, denom df = 11, p-value = 0.521
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.4282685 5.1677383
## sample estimates:
## ratio of variances 
##           1.487676

t.test(t[,1], t[,2])

## 
##  Welch Two Sample t-test
## 
## data:  t[, 1] and t[, 2]
## t = -1.7045, df = 21.186, p-value = 0.1029
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.27797435  0.02748365
## sample estimates:
## mean of x mean of y 
##  1.770459  1.895704

df <- data.frame(meta=c(rep("WT", nrow(t)), rep("KO", nrow(t))), ch=c(t[,1],t[,2]), covar=covar)
res <- ANCOVA(df)

## linear fit of dependent variable with covariates

## 
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.6653 -0.9232 -0.0703  0.7130  4.0752 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  28.4603     3.4016   8.367 2.78e-08 ***
## df[[2]]      -0.6671     1.8464  -0.361    0.721    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.659 on 22 degrees of freedom
## Multiple R-squared:  0.005897,   Adjusted R-squared:  -0.03929 
## F-statistic: 0.1305 on 1 and 22 DF,  p-value: 0.7213
## 
## y = 28 + -0.67x, r^2 = 0.0059

## 
## _____________________
## Some stats on the variables
## DV - IV:
## df[[1]]: KO
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  12.00000000   0.00000000   0.00000000   1.56026454   2.16663334   0.60636880 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##  22.74845111   1.90175955   1.89570426   0.04658797   0.10253943   0.02604527 
##      std.dev     coef.var 
##   0.16138547   0.08513220 
## ------------------------------------------------------------ 
## df[[1]]: WT
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  12.00000000   0.00000000   0.00000000   1.48426845   2.15402819   0.66975974 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##  21.24550693   1.76912376   1.77045891   0.05682350   0.12506767   0.03874692 
##      std.dev     coef.var 
##   0.19684237   0.11118155 
## Covar - IV:
## df[[1]]: KO
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  12.00000000   0.00000000   0.00000000  24.70000000  28.00000000   3.30000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 320.70000000  26.80000000  26.72500000   0.28394542   0.62495965   0.96750000 
##      std.dev     coef.var 
##   0.98361578   0.03680508 
## ------------------------------------------------------------ 
## df[[1]]: WT
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  12.00000000   0.00000000   0.00000000  24.70000000  31.50000000   6.80000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 333.00000000  27.90000000  27.75000000   0.57689083   1.26972816   3.99363636 
##      std.dev     coef.var 
##   1.99840846   0.07201472 
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  1.2641  0.273
##       22               
## _____________________
## ANOVA of DV with IV
##             Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]]      1 0.0941 0.09412   2.905  0.102
## Residuals   22 0.7127 0.03240               
## _____________________
## ANOVA of covariate with IV, independent of DV
##             Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]]      1   6.30   6.304   2.541  0.125
## Residuals   22  54.57   2.481               
## _____________________
## ANOVA of DV with IV, normalized to covariate
##             Df    Sum Sq   Mean Sq F value Pr(>F)  
## df[[1]]      1 0.0002781 2.781e-04   5.453 0.0291 *
## Residuals   22 0.0011218 5.099e-05                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## _____________________
## ANCOVA
##             Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]]      1 0.0941 0.09412   2.777  0.110
## df[[3]]      1 0.0010 0.00099   0.029  0.866
## Residuals   21 0.7117 0.03389               
## Anova Table (Type III tests)
## 
## Response: df[[2]]
##              Sum Sq Df F value  Pr(>F)  
## (Intercept) 0.24112  1  7.1143 0.01442 *
## df[[1]]     0.09035  1  2.6657 0.11743  
## df[[3]]     0.00099  1  0.0291 0.86615  
## Residuals   0.71173 21                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

check_ANCOVA(res)

RQ

Hide

Select a tab above to view

Total cycle

total.RQ <- total %>% select(contains("RQ"))
total.RQ.WT <- total.RQ %>% select_WT()
total.RQ.KO <- total.RQ %>% select_KO()

rcf.WT <- lapply(total.RQ.WT, function(x) rcf(x, "RQ", "WT", 0.6, 1.1, 0.05))
rcf.KO <- lapply(total.RQ.KO, function(x) rcf(x, "RQ", "KO", 0.6, 1.1, 0.05))

data <- bind_rows(bind_rows(rcf.WT), bind_rows(rcf.KO)) %>% split(.$id)

Models <- lapply(data, function(x){
  nplr(x$breaks, x$RQ, silent = TRUE, useLog = F)
})

overlay(Models, xlab = expression(VO[2]~(mL/min)), ylab = "relative cumulative frequency", main="RCF of WT vs KO mice RQ in complete cycle", cex.main=1.5, lwd = 3, Cols = hue_pal()(2))

print(paste0("WT EC50: ", getEstimates(Models$WT, 0.5)$x))

## [1] "WT EC50: 0.837304959775152"

print(paste0("KO EC50: ", getEstimates(Models$KO, 0.5)$x))

## [1] "KO EC50: 0.847557697979224"

t <- list(rcf.WT, rcf.KO) %>% sapply(EC50) %>% set_rownames(NULL) %>% set_colnames(c("WT", "KO"))
var.test(t[,1], t[,2])

## 
##  F test to compare two variances
## 
## data:  t[, 1] and t[, 2]
## F = 2.0112, num df = 11, denom df = 11, p-value = 0.262
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.5789754 6.9862557
## sample estimates:
## ratio of variances 
##           2.011186

t.test(t[,1], t[,2])

## 
##  Welch Two Sample t-test
## 
## data:  t[, 1] and t[, 2]
## t = -1.6606, df = 19.771, p-value = 0.1126
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.035529592  0.004046878
## sample estimates:
## mean of x mean of y 
## 0.8336410 0.8493824

df <- data.frame(meta=c(rep("WT", nrow(t)), rep("KO", nrow(t))), ch=c(t[,1],t[,2]), covar=covar)
res <- ANCOVA(df)

## linear fit of dependent variable with covariates

## 
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.5901 -0.9979  0.0213  0.8451  4.2127 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   20.790     12.043   1.726   0.0983 .
## df[[2]]        7.662     14.305   0.536   0.5976  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.653 on 22 degrees of freedom
## Multiple R-squared:  0.01287,    Adjusted R-squared:  -0.032 
## F-statistic: 0.2869 on 1 and 22 DF,  p-value: 0.5976
## 
## y = 21 + 7.7x, r^2 = 0.0129

## 
## _____________________
## Some stats on the variables
## DV - IV:
## df[[1]]: KO
##      nbr.val     nbr.null       nbr.na          min          max        range 
## 1.200000e+01 0.000000e+00 0.000000e+00 8.225666e-01 8.818145e-01 5.924789e-02 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 1.019259e+01 8.507364e-01 8.493824e-01 5.462712e-03 1.202335e-02 3.580946e-04 
##      std.dev     coef.var 
## 1.892339e-02 2.227900e-02 
## ------------------------------------------------------------ 
## df[[1]]: WT
##      nbr.val     nbr.null       nbr.na          min          max        range 
## 12.000000000  0.000000000  0.000000000  0.787188268  0.870610070  0.083421802 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 10.003692068  0.838138435  0.833641006  0.007747015  0.017051066  0.000720195 
##      std.dev     coef.var 
##  0.026836449  0.032191853 
## Covar - IV:
## df[[1]]: KO
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  12.00000000   0.00000000   0.00000000  24.70000000  28.00000000   3.30000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 320.70000000  26.80000000  26.72500000   0.28394542   0.62495965   0.96750000 
##      std.dev     coef.var 
##   0.98361578   0.03680508 
## ------------------------------------------------------------ 
## df[[1]]: WT
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  12.00000000   0.00000000   0.00000000  24.70000000  31.50000000   6.80000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 333.00000000  27.90000000  27.75000000   0.57689083   1.26972816   3.99363636 
##      std.dev     coef.var 
##   1.99840846   0.07201472 
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1   1.078 0.3104
##       22               
## _____________________
## ANOVA of DV with IV
##             Df   Sum Sq   Mean Sq F value Pr(>F)
## df[[1]]      1 0.001487 0.0014867   2.758  0.111
## Residuals   22 0.011861 0.0005391               
## _____________________
## ANOVA of covariate with IV, independent of DV
##             Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]]      1   6.30   6.304   2.541  0.125
## Residuals   22  54.57   2.481               
## _____________________
## ANOVA of DV with IV, normalized to covariate
##             Df    Sum Sq   Mean Sq F value Pr(>F)  
## df[[1]]      1 1.656e-05 1.656e-05     5.3 0.0312 *
## Residuals   22 6.874e-05 3.124e-06                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## _____________________
## ANCOVA
##             Df   Sum Sq   Mean Sq F value Pr(>F)
## df[[1]]      1 0.001487 0.0014867   2.804  0.109
## df[[3]]      1 0.000726 0.0007262   1.370  0.255
## Residuals   21 0.011135 0.0005302               
## Anova Table (Type III tests)
## 
## Response: df[[2]]
##               Sum Sq Df F value    Pr(>F)    
## (Intercept) 0.042923  1 80.9508 1.192e-08 ***
## df[[1]]     0.002041  1  3.8496   0.06315 .  
## df[[3]]     0.000726  1  1.3697   0.25498    
## Residuals   0.011135 21                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

check_ANCOVA(res)

Light cycle

light.RQ <- light %>% select(contains("RQ"))
light.RQ.WT <- light.RQ %>% select_WT()
light.RQ.KO <- light.RQ %>% select_KO()

rcf.WT <- lapply(light.RQ.WT, function(x) rcf(x, "RQ", "WT", 0.6, 1.1, 0.05))
rcf.KO <- lapply(light.RQ.KO, function(x) rcf(x, "RQ", "KO", 0.6, 1.1, 0.05))

data <- bind_rows(bind_rows(rcf.WT), bind_rows(rcf.KO)) %>% split(.$id)

Models <- lapply(data, function(x){
  nplr(x$breaks, x$RQ, silent = TRUE, useLog = F)
})

overlay(Models, xlab = expression(VO[2]~(mL/min)), ylab = "relative cumulative frequency", main="RCF of WT vs KO mice RQ in light cycle", cex.main=1.5, lwd = 3, Cols = hue_pal()(2))

print(paste0("WT EC50: ", getEstimates(Models$WT, 0.5)$x))

## [1] "WT EC50: 0.784812573550401"

print(paste0("KO EC50: ", getEstimates(Models$KO, 0.5)$x))

## [1] "KO EC50: 0.792340459568411"

t <- list(rcf.WT, rcf.KO) %>% sapply(EC50) %>% set_rownames(NULL) %>% set_colnames(c("WT", "KO"))
var.test(t[,1], t[,2])

## 
##  F test to compare two variances
## 
## data:  t[, 1] and t[, 2]
## F = 2.8464, num df = 11, denom df = 11, p-value = 0.09688
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.8194132 9.8875191
## sample estimates:
## ratio of variances 
##           2.846395

t.test(t[,1], t[,2])

## 
##  Welch Two Sample t-test
## 
## data:  t[, 1] and t[, 2]
## t = -0.56504, df = 17.88, p-value = 0.5791
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.02915937  0.01680362
## sample estimates:
## mean of x mean of y 
## 0.7858204 0.7919983

df <- data.frame(meta=c(rep("WT", nrow(t)), rep("KO", nrow(t))), ch=c(t[,1],t[,2]), covar=covar)
res <- ANCOVA(df)

## linear fit of dependent variable with covariates

## 
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.3998 -1.1543 -0.2562  1.0993  3.1561 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)    8.616      9.587   0.899   0.3786  
## df[[2]]       23.604     12.146   1.943   0.0649 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.537 on 22 degrees of freedom
## Multiple R-squared:  0.1465, Adjusted R-squared:  0.1077 
## F-statistic: 3.777 on 1 and 22 DF,  p-value: 0.06487
## 
## y = 8.6 + 24x, r^2 = 0.147

## 
## _____________________
## Some stats on the variables
## DV - IV:
## df[[1]]: KO
##      nbr.val     nbr.null       nbr.na          min          max        range 
## 1.200000e+01 0.000000e+00 0.000000e+00 7.573150e-01 8.328037e-01 7.548870e-02 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 9.503979e+00 7.943118e-01 7.919983e-01 5.574838e-03 1.227013e-02 3.729458e-04 
##      std.dev     coef.var 
## 1.931180e-02 2.438364e-02 
## ------------------------------------------------------------ 
## df[[1]]: WT
##      nbr.val     nbr.null       nbr.na          min          max        range 
## 12.000000000  0.000000000  0.000000000  0.736403245  0.843938825  0.107535580 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
##  9.429844919  0.783804817  0.785820410  0.009405455  0.020701266  0.001061551 
##      std.dev     coef.var 
##  0.032581451  0.041461701 
## Covar - IV:
## df[[1]]: KO
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  12.00000000   0.00000000   0.00000000  24.70000000  28.00000000   3.30000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 320.70000000  26.80000000  26.72500000   0.28394542   0.62495965   0.96750000 
##      std.dev     coef.var 
##   0.98361578   0.03680508 
## ------------------------------------------------------------ 
## df[[1]]: WT
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  12.00000000   0.00000000   0.00000000  24.70000000  31.50000000   6.80000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 333.00000000  27.90000000  27.75000000   0.57689083   1.26972816   3.99363636 
##      std.dev     coef.var 
##   1.99840846   0.07201472 
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value  Pr(>F)  
## group  1  3.0908 0.09264 .
##       22                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## _____________________
## ANOVA of DV with IV
##             Df   Sum Sq   Mean Sq F value Pr(>F)
## df[[1]]      1 0.000229 0.0002290   0.319  0.578
## Residuals   22 0.015779 0.0007172               
## _____________________
## ANOVA of covariate with IV, independent of DV
##             Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]]      1   6.30   6.304   2.541  0.125
## Residuals   22  54.57   2.481               
## _____________________
## ANOVA of DV with IV, normalized to covariate
##             Df    Sum Sq   Mean Sq F value Pr(>F)  
## df[[1]]      1 9.590e-06 9.594e-06   4.187 0.0529 .
## Residuals   22 5.041e-05 2.291e-06                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## _____________________
## ANCOVA
##             Df   Sum Sq  Mean Sq F value Pr(>F)  
## df[[1]]      1 0.000229 0.000229   0.381  0.544  
## df[[3]]      1 0.003169 0.003169   5.277  0.032 *
## Residuals   21 0.012610 0.000600                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
## 
## Response: df[[2]]
##                Sum Sq Df F value    Pr(>F)    
## (Intercept) 0.0262812  1 43.7657 1.501e-06 ***
## df[[1]]     0.0010525  1  1.7528   0.19976    
## df[[3]]     0.0031690  1  5.2773   0.03199 *  
## Residuals   0.0126104 21                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

check_ANCOVA(res)

## [1] "Varies significantly with Covar, p = 0.03198616938202"

Dark cycle

dark.RQ <- dark %>% select(contains("RQ"))
dark.RQ.WT <- dark.RQ %>% select_WT()
dark.RQ.KO <- dark.RQ %>% select_KO()

rcf.WT <- lapply(dark.RQ.WT, function(x) rcf(x, "RQ", "WT", 0.6, 1.1, 0.05))
rcf.KO <- lapply(dark.RQ.KO, function(x) rcf(x, "RQ", "KO", 0.6, 1.1, 0.05))

data <- bind_rows(bind_rows(rcf.WT), bind_rows(rcf.KO)) %>% split(.$id)

Models <- lapply(data, function(x){
  nplr(x$breaks, x$RQ, silent = TRUE, useLog = F)
})

overlay(Models, xlab = expression(VO[2]~(mL/min)), ylab = "relative cumulative frequency", main="RCF of WT vs KO mice RQ in dark cycle", cex.main=1.5, lwd = 3, Cols = hue_pal()(2))

print(paste0("WT EC50: ", getEstimates(Models$WT, 0.5)$x))

## [1] "WT EC50: 0.878325786494703"

print(paste0("KO EC50: ", getEstimates(Models$KO, 0.5)$x))

## [1] "KO EC50: 0.903008775451154"

t <- list(rcf.WT, rcf.KO) %>% sapply(EC50) %>% set_rownames(NULL) %>% set_colnames(c("WT", "KO"))
var.test(t[,1], t[,2])

## 
##  F test to compare two variances
## 
## data:  t[, 1] and t[, 2]
## F = 1.2216, num df = 11, denom df = 11, p-value = 0.7458
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.3516579 4.2433103
## sample estimates:
## ratio of variances 
##           1.221554

t.test(t[,1], t[,2])

## 
##  Welch Two Sample t-test
## 
## data:  t[, 1] and t[, 2]
## t = -2.5364, df = 21.783, p-value = 0.01889
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.048071153 -0.004809213
## sample estimates:
## mean of x mean of y 
## 0.8730539 0.8994941

df <- data.frame(meta=c(rep("WT", nrow(t)), rep("KO", nrow(t))), ch=c(t[,1],t[,2]), covar=covar)
res <- ANCOVA(df)

## linear fit of dependent variable with covariates

## 
## Call:
## lm(formula = df[[3]] ~ df[[2]], data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.4826 -1.1494 -0.1335  0.8566  3.8316 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    40.25      10.47   3.844 0.000882 ***
## df[[2]]       -14.68      11.81  -1.243 0.226834    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.608 on 22 degrees of freedom
## Multiple R-squared:  0.06565,    Adjusted R-squared:  0.02318 
## F-statistic: 1.546 on 1 and 22 DF,  p-value: 0.2268
## 
## y = 40 + -15x, r^2 = 0.0657

## 
## _____________________
## Some stats on the variables
## DV - IV:
## df[[1]]: KO
##      nbr.val     nbr.null       nbr.na          min          max        range 
## 1.200000e+01 0.000000e+00 0.000000e+00 8.553261e-01 9.413496e-01 8.602341e-02 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 1.079393e+01 9.058534e-01 8.994941e-01 6.993827e-03 1.539331e-02 5.869634e-04 
##      std.dev     coef.var 
## 2.422733e-02 2.693439e-02 
## ------------------------------------------------------------ 
## df[[1]]: WT
##      nbr.val     nbr.null       nbr.na          min          max        range 
## 1.200000e+01 0.000000e+00 0.000000e+00 8.269545e-01 8.981003e-01 7.114580e-02 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 1.047665e+01 8.850517e-01 8.730539e-01 7.729852e-03 1.701329e-02 7.170074e-04 
##      std.dev     coef.var 
## 2.677699e-02 3.067049e-02 
## Covar - IV:
## df[[1]]: KO
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  12.00000000   0.00000000   0.00000000  24.70000000  28.00000000   3.30000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 320.70000000  26.80000000  26.72500000   0.28394542   0.62495965   0.96750000 
##      std.dev     coef.var 
##   0.98361578   0.03680508 
## ------------------------------------------------------------ 
## df[[1]]: WT
##      nbr.val     nbr.null       nbr.na          min          max        range 
##  12.00000000   0.00000000   0.00000000  24.70000000  31.50000000   6.80000000 
##          sum       median         mean      SE.mean CI.mean.0.95          var 
## 333.00000000  27.90000000  27.75000000   0.57689083   1.26972816   3.99363636 
##      std.dev     coef.var 
##   1.99840846   0.07201472 
## levene's test
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  0.6755 0.4199
##       22               
## _____________________
## ANOVA of DV with IV
##             Df   Sum Sq  Mean Sq F value Pr(>F)  
## df[[1]]      1 0.004194 0.004194   6.433 0.0188 *
## Residuals   22 0.014344 0.000652                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## _____________________
## ANOVA of covariate with IV, independent of DV
##             Df Sum Sq Mean Sq F value Pr(>F)
## df[[1]]      1   6.30   6.304   2.541  0.125
## Residuals   22  54.57   2.481               
## _____________________
## ANOVA of DV with IV, normalized to covariate
##             Df    Sum Sq   Mean Sq F value Pr(>F)  
## df[[1]]      1 2.661e-05 2.661e-05   5.489 0.0286 *
## Residuals   22 1.066e-04 4.847e-06                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## _____________________
## ANCOVA
##             Df   Sum Sq  Mean Sq F value Pr(>F)  
## df[[1]]      1 0.004194 0.004194   6.237 0.0209 *
## df[[3]]      1 0.000220 0.000220   0.327 0.5734  
## Residuals   21 0.014124 0.000673                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Anova Table (Type III tests)
## 
## Response: df[[2]]
##               Sum Sq Df  F value    Pr(>F)    
## (Intercept) 0.068979  1 102.5628 1.552e-09 ***
## df[[1]]     0.003197  1   4.7543   0.04075 *  
## df[[3]]     0.000220  1   0.3272   0.57337    
## Residuals   0.014124 21                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

check_ANCOVA(res)

## [1] "Varies significantly with DV, p = 0.0407456652336419"

PCA

Libraries

reqpkg <- c("ggfortify","rgl","lfda")
for (i in reqpkg) {
  print(i)
  print(packageVersion(i))
  suppressWarnings(suppressMessages(library(i, quietly=T, verbose=T, warn.conflicts=F,character.only=T)))
}

## [1] "ggfortify"
## [1] '0.4.10'
## [1] "rgl"
## [1] '0.100.54'
## [1] "lfda"
## [1] '1.1.3'

knitr::knit_hooks$set(webgl = hook_webgl)

Functions

pca_mutate <- function(data, channel) {
  df <- data %>% select(contains(channel))
  WT <- df %>% select_WT() %>% set_colnames(paste0("WT_", colnames(.))) %>% melt() %>% select(-variable)
  KO <- df %>% select_KO() %>% set_colnames(paste0("KO_", colnames(.))) %>% melt() %>% select(-variable)
  df <- bind_cols(WT, KO) %>% set_colnames(c("WT","KO")) %>% melt()
  
  return(df)
}

plot_pca <- function(pca, data, title="") {
  autoplot(pca, data = na.omit(data), colour = "id", loadings = TRUE, loadings.colour = "black", loadings.label = TRUE, loadings.label.colour = "black", loadings.label.size = 5) + ggtitle(title)
}

Light

Small filtering step: during the day, there were points where the gas analyzers failed, leaving VO2 and RQ values unbelievably low. The rows contianing data for those points were simply discarded.

light <- light %>% slice(-320:-325)

VO2_light <- pca_mutate(light, "VO2")
RQ_light <- pca_mutate(light, "RQ")
Food_light <- pca_mutate(light, "Food")
Water_light <- pca_mutate(light, "Water")
BM_light <- pca_mutate(light, "BodyMass")
PedMeters_light <- pca_mutate(light, "PedMeters_R")

df_light <- RQ_light %>% set_colnames(c("id","RQ")) %>% add_column(VO2=VO2_light$value) %>% add_column(FoodIn=Food_light$value) %>% add_column(WaterIn=Water_light$value) %>% add_column(BodyMass=BM_light$value) %>% add_column(PedMeters=PedMeters_light$value)

PCA of all variables

df_light %>% select(-id) %>% prcomp(scale. = TRUE) %>% plot_pca(df_light, title = "PCA of all variables, light cycle")

PCA without BodyMass

df_light %>% select(-id,-BodyMass) %>% prcomp(scale. = TRUE) %>% plot_pca(df_light, title="PCA of all variables besides BM, light cycle")

VO2 vs RQ

df_light %>% select(VO2,RQ) %>% prcomp(scale. = TRUE) %>% plot_pca(df_light, title="PCA VO2 and RQ, light cycle")

Food vs Water

df_light %>% select(FoodIn,WaterIn) %>% prcomp(scale. = TRUE) %>% plot_pca(df_light, title="PCA Food and Water, light cycle")

Dark

VO2_dark <- pca_mutate(dark, "VO2")
RQ_dark <- pca_mutate(dark, "RQ")
Food_dark <- pca_mutate(dark, "Food")
Water_dark <- pca_mutate(dark, "Water")
BM_dark <- pca_mutate(dark, "BodyMass")
PedMeters_dark <- pca_mutate(dark, "PedMeters_R")

df_dark <- RQ_dark %>% set_colnames(c("id","RQ")) %>% add_column(VO2=VO2_dark$value) %>% add_column(FoodIn=Food_dark$value) %>% add_column(WaterIn=Water_dark$value) %>% add_column(BodyMass=BM_dark$value) %>% add_column(PedMeters=PedMeters_dark$value)

PCA of all variables

df_dark %>% select(-id) %>% prcomp(scale. = TRUE) %>% plot_pca(df_dark, "PCA of all variables, dark cycle")

PCA without BodyMass

df_dark %>% select(-id,-BodyMass) %>% prcomp(scale. = TRUE) %>% plot_pca(df_dark, title="PCA of all variables besides BM, dark cycle")

VO2 vs RQ

df_dark %>% select(VO2,RQ) %>% prcomp(scale. = TRUE) %>% plot_pca(df_dark, title="PCA VO2 and RQ, dark cycle")

Food vs Water

df_dark %>% select(FoodIn,WaterIn) %>% prcomp(scale. = TRUE) %>% plot_pca(df_dark, title="PCA Food and Water, dark cycle")

Local Fisher Discriminant Analysis (LDFA)

I attempted to model and plot ldfa. LDFA should give similar results to PCA, with the advantage of being sensitive to multimodality of data, i.e. multiple effects driving an observation, as well as sensitive to distinctiveness in the data that might be lost in reducing dimensions in PCA. In other words, a learning model like LDFA might account for more relationships than a principal component.

More information:

Tang Y, Cohen J, Local Fisher Discriminant Analysis on Beer Style Clustering.
Yang Y, Li W, lfda: An R Package for Local Fisher Discriminant Analysis and Visualization. The R Journal Vol. XX/YY, AAAA 20ZZ

Light

model <- lfda(df_light[-1], df_light[, 1], r=3, metric = "plain")
autoplot(model, data=df_light, frame = T, colour = "id")

plot(model, labels = df_light[, 1])

You must enable Javascript to view this page properly.

Dark

model <- lfda(df_dark[-1], df_dark[, 1], r=3, metric = "plain")
autoplot(model, data=df_dark, frame = T, colour = "id")

plot(model, labels = df_dark[, 1])

You must enable Javascript to view this page properly.

Session info

## R version 3.6.3 (2020-02-29)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 18362)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=English_United States.1252 
## [2] LC_CTYPE=English_United States.1252   
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.1252    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] lfda_1.1.3       rgl_0.100.54     ggfortify_0.4.10 WRS2_1.0-0      
##  [5] pastecs_1.3.21   multcomp_1.4-12  TH.data_1.0-10   MASS_7.3-51.5   
##  [9] survival_3.1-11  mvtnorm_1.1-0    effects_4.1-4    compute.es_0.2-4
## [13] car_3.0-7        carData_3.0-3    nplr_0.1-7       purrr_0.3.3     
## [17] scales_1.1.0     tibble_2.1.3     docstring_1.0.0  reshape2_1.4.3  
## [21] DT_0.13          ggpubr_0.2.5     ggplot2_3.3.0    dplyr_0.8.5     
## [25] magrittr_1.5    
## 
## loaded via a namespace (and not attached):
##  [1] minqa_1.2.4             colorspace_1.4-1        ggsignif_0.6.0         
##  [4] rio_0.5.16              htmlTable_1.13.3        markdown_1.1           
##  [7] base64enc_0.1-3         mc2d_0.1-18             rstudioapi_0.11        
## [10] roxygen2_7.1.0          farver_2.0.3            RSpectra_0.16-0        
## [13] xml2_1.2.5              codetools_0.2-16        splines_3.6.3          
## [16] knitr_1.28              Formula_1.2-3           jsonlite_1.6.1         
## [19] nloptr_1.2.2.1          cluster_2.1.0           png_0.1-7              
## [22] shiny_1.4.0.2           compiler_3.6.3          backports_1.1.5        
## [25] fastmap_1.0.1           assertthat_0.2.1        Matrix_1.2-18          
## [28] survey_3.37             later_1.0.0             acepack_1.4.1          
## [31] htmltools_0.4.0         tools_3.6.3             gtable_0.3.0           
## [34] glue_1.3.2              Rcpp_1.0.4.6            cellranger_1.1.0       
## [37] vctrs_0.2.4             nlme_3.1-144            crosstalk_1.1.0.1      
## [40] xfun_0.12               stringr_1.4.0           openxlsx_4.1.4         
## [43] lme4_1.1-21             miniUI_0.1.1.1          mime_0.9               
## [46] lifecycle_0.2.0         klippy_0.0.0.9500       zoo_1.8-7              
## [49] hms_0.5.3               promises_1.1.0          sandwich_2.5-1         
## [52] RColorBrewer_1.1-2      yaml_2.2.1              curl_4.3               
## [55] gridExtra_2.3           rpart_4.1-15            reshape_0.8.8          
## [58] latticeExtra_0.6-29     stringi_1.4.6           checkmate_2.0.0        
## [61] manipulateWidget_0.10.1 boot_1.3-24             zip_2.0.4              
## [64] rlang_0.4.5             pkgconfig_2.0.3         evaluate_0.14          
## [67] lattice_0.20-40         htmlwidgets_1.5.1       labeling_0.3           
## [70] cowplot_1.0.0           tidyselect_1.0.0        plyr_1.8.6             
## [73] R6_2.4.1                Hmisc_4.4-0             DBI_1.1.0              
## [76] pillar_1.4.3            haven_2.2.0             foreign_0.8-75         
## [79] withr_2.1.2             mgcv_1.8-31             abind_1.4-5            
## [82] nnet_7.3-13             crayon_1.3.4            rARPACK_0.11-0         
## [85] rmarkdown_2.1           jpeg_0.1-8.1            grid_3.6.3             
## [88] readxl_1.3.1            data.table_1.12.8       forcats_0.5.0          
## [91] webshot_0.5.2           digest_0.6.25           xtable_1.8-4           
## [94] tidyr_1.0.2             httpuv_1.5.2            munsell_0.5.0          
## [97] mitools_2.4

By Sumeed Yoyo Manzoor

smanzoor@uchicago.edu

RCF and PCA of WT vs L-Bmal1-KO RMR data from metabolic cages

Sumeed ManzoorChang Lab at the University of Chicago

May 23, 2020, updated May 28, 2020

Dataset

Create data

Import data

Datatables

Light

Hide

VO2

RQ

BodyMass

Food

Water

PedMeters

Dark

Hide

VO2

RQ

BodyMass

Food

Water

PedMeters

Graphs

Total cycle

Hide

VO2

RQ

BodyMass

Food

Water

PedMeters

RCF

Libraries

Functions

VO2

Hide

Total cycle

Light cycle

Dark cycle

RQ

Hide

Total cycle

Light cycle

Dark cycle

PCA

Libraries

Functions

Light

Dark

Local Fisher Discriminant Analysis (LDFA)

Light

Dark

Session info

Sumeed Manzoor
Chang Lab at the University of Chicago